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1. INTRODUCTION 

In this paper we shall describe a new numerical method for solving conservation laws. It is 
much simpler than a typical high resolution method [1]. No flux limiter or any characteristics- 
based technique is involved. No artificial viscosity or smoothing is introduced, and no moving 
mesh is used. Yet this method is capable of generating highly accurate shock tube solutions. The 
slight numerical overshoot and/or oscillations generated can be removed if a simple averaging 
formula initially used is replaced by a weighted averaging formula. This modification has no 
discemable effect on other parts of the solution. Because of its simplicity, multi-dimension 
generalization is straightforward and it allows for the simultaneous treatment of variables in 
different spatial directions. 


2. CONSERVATION LAWS 


We consider a dimensionless form of the 1-D unsteady Euler equations for an ideal gas. Let 
p, u, p, and y, respectively, be the mass density, velocity, static pressure, and constant specific 
heat ratio. Let 


q\ =p , <72 = P“ < < ?3 =p/(Y-l) + (l/2)p« 2 

/l=<?2 

fi = (Y- 1) <73 + (1/2) (3-y) (<? 2 ) 2 /<7 l 
h = Y<?2<?3/<7i ~ (1/2) (Y~ 1) (<7 2> 3 /(<? l ) 2 . 

Then the Euler equations can be expressed as 

dq m /dt+df m /dx = 0 , m = 1, 2, 3 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 


§ This work is dedicated to the memory of a teacher, Mr. Nylon Cheng 
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Let x i = x and x 2 = r be considered as the coordinates of a two-dimensional Euclidean space 
E 2 . The integral form of Eq. (2.3) in the space-time £ 2 can be expressed as (see Fig. 1) 

<6 S..2-0 , m = 1, 2, 3 (2.4) 

7 S(V) m 

where (i) S(V) is the boundary of an arbitrary space-time volume V in E 2 , (ii) X* = (/*.?«.) are 
space-time current density vectors in E 2 , and (iii) ds — dolt with do and it, respectively, being 
the area and the outward unit normal of a surface element on S(V). Note that (i) h m -ds is the 
space-time flux of X m leaving the volume V through the surface element ds, and (ii) all 
mathematical operations can be carried out as though E 2 is an ordinary two-dimensional 
Euclidean space. 


3. NUMERICAL METHOD 

Let E 2 be divided into nonoverlapping rhombic regions (see Fig. 2) referred to as solution 
elements (SEs). Each SE is centered at a mesh point (j,n) where n = 0, 1/2, 1, 3/2, • • • , and j = 
(/i ±1/2), (n±3/2), • • • , i.e., j is a half-integer (whole integer) if n is a whole integer (half-integer). 
In other words, j and n are both whole integers or both half-integers if, as occurring in Fig. 2, a 
SE is centered at the mesh point (j,n+ 1/2). Thus, the locations of SEs and their centers are 
staggered over every half time-step. A SE centered at (J,n), and its interior are denoted by 
SE (j,n) and SE' (J,n), respectively. 

For any (x,t) eSE'O'.n), q m (x,t), f m (x,t), and Ji m (x.t), respectively, are approximated by 
qjjc, r,j,n),f m (x,t;j,n), and Xn(*>/ ;/',«) which we shall define immediately. Let 

q m (x,r,j,n) = (o m )" + (a*)" (x -Xj) + (p m )y (t-t n ) , m = l,2, 3 (3.1) 

where (a m )", (o^)", and (p OT )" are constants in SE'(J,n), and (xj,t n ) are the coordinates of the 
mesh point (j,n). Note that 

q_ m (Xj,t n J,n) = (o m )j , dq m (x,rj,n)/dx = (a m )] , dq m (x,c,j,n)/dt = (p m )" (3.2) 

Moreover, if we identify (o m )", (o^,)", and (P m );, respectively, with the values of q m , dq m /dx, 
and dq m /Bt at ( xj,t n ), the expression on the right side of Eq. (3.1) becomes the first-order 
Taylor’s expansion of q m (x,t ) at (x j,t n ). As a result of these considerations, (o m )", (a m )", and 
(P m )" will be considered as the numerical analogues of the values of q m , dq m /dx, and dq m /dt at 
( xj,t n ), respectively. 

Let f m (x,t;j,n), m= 1, 2, 3, be defined in terms of q m (x,f,j,n), m= 1, 2, 3, according to Eq. 
(2.2) with the understanding that f m and q m in Eq. (2.2) be replaced by f m (x,t;j,n ) and 
q_ m (x,t\j,n), respectively. Using Eq. (3.1), f m (x,t\j,n ) is expressed as a function of (x-Xj) and 
(r-r*), and then expanded as a power series of them. The new method is simplified by truncating 
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the power series after first-order terms. This is consistent with the first-order approximation 
given in Eq. (3.1). In the new method, only f m at x = x, arc needed. Let o m = (a m )" and p m = 
(p m )". Then explicitly, we have: 

fi(Xj,r,j,n) = c 2 + $ 2 (t-t n ) (3.3) 

) = (y- 1 )o 3 + (1/2)(3 -y)(o 2 ) 2 /o, 

+ [ (Y“ 1) p3 + (3 ■ -Y) (O 2 P 2 /Oi) - (1/2) (3-Y)(o 2 /0i) 2 Pi ] (t-t n ) (3.4) 

/ 3 (x ; ,r;7,n) = Y0 2 o 3 /o 1 -(1/2)(y-1)(o 2 ) 3 /(0!) 2 + {YK 02 P 3 +<J3p2)/Oi 

- a 2 o 3 pi /(o,) 2 ] + (l/2)(Y-l)[2(a 2 /o,) 3 pi -3(o 2 /o,) 2 p 2 ] } (t-t n ) (3.5) 
Since X. = we define ~i m i.x,t\j,n) = (f m (x,t\j,n ) , q m (x,fj,n)). 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 2) referred to as 
conservation elements (CEs). They are also staggered over every half time-step. Let the CE with 
its upper edge centered at ( j,n ) be denoted by CE (J,n). Then the current approximation of Eq. 
(2.4) is 

<6 lL-ds= 0 (all possible m and (J,n) ) (3.6) 

"S(CE0») 

Because the entire boundary (except for three isolated points) of a CE is located within the 
interiors of three neighboring SEs, X is continuous across any interface separating two 
neighboring CEs. Thus Eq. (3.6) will remain valid if CE(y,n) is replaced by the union of any 
combination of CEs. 

Because each S( CEO', «)) is a simple closed curve in £ 2 (see Fig. 1), the surface integration 
form Eq. (3.6) can be converted into a line integration form [2, p.14], i.e., 

<6 X • dr = 0 (all possible m and (j,n) ) (3.7) 

7s(ceo.«)) ~ 

where X = (- £m - fm ) ar >d dr = ( dx,dt ). 

For each SE(J,n), let 

s,0',n) = (Ax/8)a, +(l/2)(Ar/Ax)[o 2 + (Ar/4)P2 ] (3.8) 

j 2 0',n) = (Ax/8)a 2 + (l/2)(Ar/Ax)((Y-l)o 3 +(l/2)(3-Y)(o 2 ) 2 /Oi ] 

+ (l/8)[(Ar) 2 /Ax] {(Y-l)p3 + (3-Y)[02p2/0t - (l/2)(a 2 /a,) 2 pi ] } (3.9) 

s 2 (j,n) = (Ax/8) a 3 + (l/2)(Ar/ Ax) [ y<J 2 o 3 /a, - (l/2)(y- l)(o 2 ) 3 /(o,) 2 ] + (1/8) [ (Ar) 2 /Ax ] 
x {(Y/°i)(° 2P3 +O3P2 -o 2 o 3 Pi/Oi) + (y- 1 )[(o 2 /Oi) 3 Pi -(3/2)(a 2 /0!) 2 P2] } (3.10) 
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where o m ~ (o m )", a m = (a*)*, and p m = (p m )". Then Eq. (3.7) implies that, for each 
SE(/,rt+l/2), 

(Om)" +1/2 = < 1/2 ) t(Om)"-i/2 + (Om)"+i/2 1 + .r,„0'-l/2,n) - j m 0+l/2,rt) , m = 1, 2, 3 (3.1 1) 

i.e., (o m )" +1/2 is determined in terms of the numerical variables associated with SE(y-l/2,n) and 
SE(/+l/2,rt). Similar formulae for (a m )" +I/2 and (p m )" +,/2 wiU be given next. 

Let i4 + . A, and A _ (see Fig. 2) denote (x, + i/ 2 ,r" +1/2 ), (x;,r" +1/2 ), and (x ; _i /2 ,r' ,+1/2 ), 
respectively. Let 

q_m± =q m (.x j±v2 ,t n+m -,j±\l2,n) (3.12) 

Because A ± do not belong to SE'(/±l/2,n), the expression on the right side of Eq. (3.12) is to be 
evaluated at the two points immediately below them. A central-difference formula for evaluating 
(a m )" +1/2 , the numerical analogue of dq m /dx at point A, is 

(a^ 172 = (q m+ - q_ m .)l Ax (3.13) 

This formula is valid as long as no discontinuity of q m (or its derivatives) occurs between A_ and 
A+. In the following discussion, we develop an alternate which is valid even in the presence of 
discontinuity. 

Let 

( 0 U±)r 1/2 = ± l^±-(o„)r 1/2 ]/(Ax/ 2 ) (3.14) 

where (a m )" +1/2 has just been obtained using Eq. (3.11). Because q m +, (a m )" +1/2 , and q_ m - are 
the numerical analogues of q m at A + , A, and A_, respectively, (a m+ )" +,/2 and are two 

numerical analogues of dq m (Xj,t n+u2 )/dx with one being evaluated from the right and another 
from the left. Note that (a m )" +1/2 defined by Eq. (3.13) is equal to the average of (a m+ )* +1/2 and 
(a m _); +1/2 . 

In case that a discontinuity occurs between A and A+ but not between A and A_, one would 
expect that |(a m+ )" +1/2 1 » |(a m _)" +1/2 |. Moreover, because A and A_ are on the same side of the 
discontinuity while A and A+ are on the opposite sides, (ot m )" +l/2 should be closer to (a m _)" +1/2 
than (a m+ )" +l/2 . This observation suggests that (a m )" +1/2 should be a weighted average of 
(a m+ )" +1/2 and (a m _)J +U2 biased toward the one with the smaller magnitude. 

As a result of the above and other considerations [3], Eq. (3.13) will be generalized by 

(a m )" +1/2 =F((a m _)" +1/2 ,(a m+ )" +I/2 ;c) (3.15) 

Here c>0 is an adjustable constant and the function F is defined by (i) F(0,0;c) = 0 and (ii) 
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F(cx_,a + ;c) = (|a + | c CL + |cc| c a f )/(|a + | c + |ol| c ) , (la, |+ |a_| > 0) (3.16) 

where 04 and 04 are any two real variables. Note that F(cl,C 4 ;c) = (cl + 04 ) /2 if c = 0 or |cl | 
= | 04 |, i.e., Eq. (3.15) is reduced to Eq. (3.13) if c = 0 or |ol| = |cs*|. Also the expression on the 
right side of Eq. (3.16) represents a weighted average of cl and 04 with the weight factors 
| 04| c /(|04 | c + 1 04 | c ) and |c 4 | c /(I 04 | c + |ol | c ). For c > 0, this average is biased toward the one 
among cl and 04 with the smaller magnitude. For the same values of [ 04 1 and |cl|, the bias 
increases as c increases. 

Substituting Eq. (2.2) into Eq. (2.3), one obtains three equations in which dq m /dt, m = 1, 2, 3, 
are expressed in terms of q m and dq m / dx, m = 1, 2, 3. Let q m and their derivatives be replaced by 
the corresponding numerical analogues at the mesh point ( j.n+ 1 / 2 ), one obtains that 

(P,)r ,/ 2 =-a 2 (3-17) 

(P2)j + 1 / 2 =(l/2)(3-Y)(o 2 /Oi) 2 a, -(3-Y)(a 2 /ai)a 2 -(Y-l)a3 (3-18) 

(fc )" +1/2 = [ YC 2 O 3 / (ai ) 2 - (y- 1 ) (o 2 /01 ) 3 ] a, 

+ [(3/2)(Y-l)(a 2 /Oi) 2 -Y0 3 /o 1 ]a 2 -Y(o 2 /o,)a 3 (3.19) 

where o m = (G m )" +1/2 and cl, = (a m )” +1/2 . 

With the aid of Eqs. (3.11), (3.15), and (3.17) - (3.19), ( o m )J , «x m )J, and (p m )" can be 
determined in terms of the initial values (o m )±i/ 2 , (o m )± 3 / 2 . ' ' ‘ » 3°^ (ctm)±i/2> ( a m)± 3 / 2 - 


4. NUMERICAL RESULTS 

We consider a shock tube problem used by Sod [4], Let Y =1-4. At t = 0, let (i) (p ,u,p) = 
(1,0,1), i.e., 0 ?i,< 7 2 ,<7 3 ) = (1, 0,2.5), if x < 0, and (ii) (p ,u,p) = (0.125,0,0.1), i.e., (<7 i,< 7 2 ,<? 3 ) = 
(0.125,0,0.25), if x > 0. Thus 


O') 


(( 0 ])°,( 0 2 )®,( 0 3 )°) 


(1, 0,2.5) if -1/2, -3/2, 

(0.125,0,0.25) if j = 1/2, 3/2, - - - 


(4.1) 


and (ii) (c O? = 0, j = ±1/2, ±3/2, • • • . Eqs. (3.17) - (3.19) imply that (p m )? = 0, ; = ±1/2, ±3/2, 


The above initial conditions, and Eqs. (3.11), (3.15), and (3.17) — (3.19), imply that (o m )", 
( 04 ,)", and (P m )" are constant in two separate regions which, respectively, are defined by j < 
-(n+1/2) and j > (n+1/2). Thus one needs to evaluate the above variables only if (n+1/2) > \j |. 

The current scheme is stable if CFL = max (|«| + |a|)Ar/Ax < 1 [3]. Here a = local sound 
speed. In the current computations, Ax = 0.01 , Ar = 0.004, and CFL = 0.88. Numerical results 
(dots) at r = 0.4 are compared with the exact solutions (solid lines) in Fig. 3. Since each 
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marching step advances the solution from t to t+ At/2, these results are obtained after 200 steps. 
Note that (i) shock discontinuity is resolved almost within one mesh interval, and (ii) the slight 
numerical overshoot and oscillations generated when c = 0 essentially disappear when c = 1 is 
used. 


5. CONCLUSIONS AND DISCUSSIONS 

The current scheme has a stencil containing only two points. This minimization of stencil has 
the effect of reducing numerical diffusion [5]. It is achieved by including (a,*)" and ((3 m )" as 
numerical variables. The fluxes at an interface separating two CEs are evaluated with no 
interpolation or extrapolation. Accuracy of flux evaluation is enhanced by requiring that the 
solution given in Eq. (3.1) satisfies the Euler equations at the center of every SE. This makes the 
use of characteristics-based techniques less necessary. The above key features all contribute to 
the simplicity, generality, and accuracy of the current scheme. They all owe their existence to the 
use of staggered SEs and CEs. 

In the current method, flux evaluation within each SE(j,n ) is required only at a subset of 
SE(j,n), i.e., a horizontal line segment centered at (j,n) and a vertical line segment starting 
upward from (j,n) (see Fig. 2 ). As a result, we may redefine SEO'.n) to be this subset. This new 
definition is used in the following sketch of an extension of the the current scheme to a three- 
dimensional Euclidean space E 3 (xj = x, *2 = >'. and *3 = r). 

A SE contains three mutually perpendicular rectangles (see Fig. 4a). The point of intersection 
is referred to as the center of this SE. The SEs are staggered in both x- and y- directions over 
every half time-step. The CEs are rectangular boxes (see Fig. 4b) also staggered in both x- and 
y- directions over every half time-step. From Fig. 4b, it is seen that the boundary of a CE can be 
divided into five parts which, respectively, belong to five neighboring SEs. As a result, the 
solution procedure described in Section 3 can be easily extended to £3. 
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Figure 1 . — A surface element ds and a line 
segment dr on the boundary S(V) of a 
volume V in a space-time E 2 . 
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Figure 2 .— The SEs (the rhombuses formed by solid lines) 
and the CEs (the rectangles formed by dashed lines). 
The dots represent the centers of SEs. 
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c 


c 


c 


c 


implicit real*8 (a-h, o-z) 


dimension 

* 

* 

* 


ql (1000) , q2 (1000) , 
q3x(1000) , qlt(1000) 
q2n(1000) , q3n(1000) 
xx(1000) 


q3 (1000) , qlx(1000) , q2x(1000) , 
q2t(1000), q3t(1000) , qln(1000) , 
, Sl(1000) f s2 (1000) , S3 (1000) , 


it = 200 
iw = 200 
iwi = 10 
dt = 0 . 4d-2 
dx = 0.1d-l 
qa = 1 . 4d0 
rhol = l.dO 
ul = O.dO 
pi = l.dO 
rhor = 0.l25d0 
ur = O.dO 
pr = 0 . Ido 
ic = 1 

open (unit=8,file='for008') 
write (8,10) it, iw, iwi, ic 
write (8,20) dt,dx,ga 
write (8,30) rhol, ul, pi 
write (8,40) rhor , ur , pr 

dtl = 0 . 5d0*dt 
dt2 = dt**2/8.d0 
hdx = 0 . 5d0*dx 
dxl = 0 . 25d0*dx 
dx2 = 0 . 125d0*dx 
gal = ga - l.dO 
ga2 = 3. do - ga 
ga3 = 0.5d0*ga2 
ga4 = 0.5d0*gal 
ql(l) = rhol 
q2(l) = rhol*ul 

q3(l) = pl/gal + 0 . 5 d 0 *rhol*ul **2 
ql (2 ) = rhor 
q2r = rhor*ur 

q3r = pr/gal + 0 . 5 d 0 *rhor*ur **2 

q2 (2) = q2r 

q3 (2) = q3r 

qlx(l) = O.dO 

qlx ( 2 ) = O.dO 

q2x(l) = O.dO 

q2x(2) = O.dO 

q3x(l) = O.dO 

q3x(2) = O.dO 

qlt(l) = O.dO 

qlt (2) = O.dO 

q2t (1) = O.dO 

q2t(2) = O.dO 

q3t ( 1) = O.dO 

q3t (2) = O.dO 

m = 2 

do 600 i = 1, it 


c 
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\ 


c 

c 


c 


c 

100 

c 


200 


do 100 j = 1 
vl = ql(j) 
v2 = q2(j) 
v3 = q3(j) 
vlt = qlt(j) 
v2t = q2t ( j ) 
v3t = q3t ( j ) 


sl(j) = dx2*qlx( j ) + (dtl*v2 + dt2*v2t)/dx 


s2 ( j ) 

* 


dx2*q2x( j ) + (dtl* (gal*v3 + (ga3/vl) *v2**2) + dt2*(gal*v3t 
+ (ga3/vl) *(2.d0*v2*v2t - v 2 ** 2 *vlt/vl) ) ) /dx 


S3(i) - dx2*q3x ( j ) + (dtl* (ga*v2*v3/vl - ga4*v2**3/vl**2) + d *2* 
t ( (ga/vl) * (v2*v3t + V3*v2t - v2*v3*vlt/vl) - (ga4*v2/vl**2) 

k * ( 3 . d 0 *v 2 *v 2 t - 2 .dO*v 2 ** 2 *vlt/vl)))/dx 


continue 


- sl(j+l) 

- s2 ( j+1) 

- S3 (j+1) 


mm = m — 1 

do 200 j = l,mn 

qln(j+l) = 0. 5d0* (ql ( j ) + ql(D+l)) + sl(]) 

q2n(j+l) = 0 . 5 d 0 *(q 2 (j) + q2(j+l)) + s2(}) 

q3n( j+1) = 0 • 5d0* (q3 ( j ) + q3(j+l)) + s3(}) 

vlxl = (qln( j+1) - ql(j) - dtl*qlt( j ) )/hdx 
vlxr = (ql ( j+1) + dtl*qlt ( j+1) - qln(;j+l) )/hdx 
v2xl = (q2n ( j+1) - q2(j) - dtl*q2t ( j) )/hdx 
v2xr = (q2 (j+1) + dtl*q2t(j+l) - q2n( j+1) )/hdx 
v3xl = (q3n(j+l) - q3(j) - dtl*q3t( j) )/hdx 
v3xr = (q3 (j+1) + dtl*q3t (j+1) - q3nO+l) )/hdx 
alx(j+l) = (vlxl* (dabs (vlxr) )**ic + vlxr* (dabs (vlxl) ) **ic)/ 

* ( (dabs (vlxl) )**ic + (dabs (vlxr) ) **ic + l.d-60) 
a2x( j+1) = (v2xl* (dabs (v2xr) ) **ic + v 2 xr*(dabs(v 2 xl) ) **xc)/ 

* ( (dabs(v2xl) )**ic + (dabs(v2xr) ) **ic + l.d-60) 
a3x ( j+1) = (v3xl* (dabs (v3xr) ) **ic + v3xr* (dabs(v3xl) )**ic)/ 

* ( (dabs(v3xl) ) **ic + (dabs(v3xr) ) **ic + l.d-60) 
continue 


c 


300 

c 


do 300 j = 2 ,m 
ql( j ) = qln(j) 
q2 ( j ) = q2n ( j ) 
q3 ( j ) = q3n(j) 
vl = ql ( j ) 
v2 = q2 ( j ) 
v3 = q3 (j) 
vlx = qlx(j) 
v2x = q2x ( j ) 
V3x = q3x( j ) 
qlt(j) = -V2x 
q2t ( j ) = 
q3t(j) = 


continue 


qa3* (v 2 /vl) **2*vlx - ga2* (v2/vl) *v2x - gal*v3x 
(ga*v2*v3/vl**2 - gal* (v2/vl) **3) *vlx + (1.5d0*gal* 
(v2/vl)**2 - ga*v3/vl) *v2x - ga* (v2/vl) *v3x 


m = m + 1 
ql(m) = rhor 
q2(m) = q2r 
q3(m) = q3r 
qlx(m) = O.dO 
q2x(m) = O.dO 
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c 


400 

c 


500 

600 

c 

10 

20 

30 

40 

50 

60 


q3x(m) = 0 . do 
qlt(m) = O.dO 
q2t(m) = O.dO 
q3t(m) = O.dO 

if (i.ne.iw) goto 600 
iw = iw + iwi 
t = dtl*dfloat (i) 
write (8/50) i,t 
mm = m — 1 
t2 = dx*dfloat (mm) 

XX (1) = -0.5d0*t2 

do 400 j = l,mm 
xx(j+l) = xx(j) + dx 
continue 


do 500 j = l,m 

X = q2(j)/ql(j) 

y = q3(j)/ql(j) - 0 . 5d0*x**2 

z = gal*y*ql(j) 

write ( 8 , 60 ) xx(j),y,ql(j),x,z 

continue 

continue 


close (unit=8) 
format ( * it = • ,U, • iw = 
format (' dt = ' ,gl4.7 , ' dx 
format (' rhol = , ,gl4.7, / 
fonnat(' rhor = , ,gl4.7, / 
format ( ' i= ' , i4 , ' t= ' , 
f ormat ( ' x =, ,f8.4, / e — , 


/ ,i4, / iwi = ' ,i4, # ic - / ,i4) 

= / ,gl4.7, / gamma = / ,gl4.7) 
ul = • pl = ',gl4.7) 
ur = ',gl4.7,' pr = ',gl4.7) 

g 14 * 7 , '*************************) 
gll.4,' rho = , ,gll.4, , u =',gll.4, 


* ' P =' ,gll.4) 


stop 

end 
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